##########################################################################################

library(ggplot2)
library(tidyverse)
library(magrittr)
library(parallel)
library(data.table)
library(rjags)
library(ggmcmc)
library(kableExtra)

##########################################################################################

input_file <- "~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime/ipmn-timing-master/data/Timing_Metrics.xlsx"
out_path <- "~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime"
code_path <- "~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime/ipmn-timing-master"

##########################################################################################

source(paste0(code_path, "/code/functions.R"))
jag_file <- paste0( code_path , "/code/clock.jag"  )

##########################################################################################


dat <- read_excel(input_file,
                sheet=2) %>%
"["(-nrow(.), ) %>%
set_colnames(unixFriendly(.)) %>%
mutate(additional=mutations_ca-mutations_hg_ipmn) %>%
filter(case != "MTP19") %>%
mutate(high_mutation_rate=additional > 20) %>%
mutate(jags_id=paste0("T[", seq_len(nrow(.)), "]"))

y <- dat$additional

##########################################################################################
## mutations/year
## 1-10

posteriors <- readRDS(paste0(code_path, "/output/timing_master.R/timing.rds")) %>%
  mutate(patient=as.character(patient)) %>%
  left_join(select(dat, jags_id, case, high_mutation_rate),
            by=c("patient"="jags_id"))

high.mr <- filter(posteriors, high_mutation_rate, mut.per.yr > 3)
low.mr  <- filter(posteriors, !high_mutation_rate, mut.per.yr < 7)

##########################################################################################

high.mr3  <- filter(posteriors, high_mutation_rate, mut.per.yr > 3)
low.mr2  <- filter(posteriors, !high_mutation_rate, mut.per.yr < 7)


high.mr3 <- smoothPosterior(high.mr)
low.mr2 <- smoothPosterior(low.mr)

#posteriors2 <- smoothPosterior(posteriors)

##########################################################################################
post <- bind_rows(low.mr2, high.mr3) %>%
  mutate(case=factor(case, levels=dat$case))
#post <- posteriors2 %>%
#  mutate(case=factor(case, levels=dat$case))

post2 <- post %>%
  group_by(case) %>%
  summarize(lower=quantile(median, 0.05),
            mean=mean(median),
            upper=quantile(median, 0.95)) %>%
  mutate(x=11)
one.number <- mean(post2$mean)
one.number <- mean(post2$mean)

p <- ggplot(post,
       aes(ymin=lower, ymax=upper,
           y=median, x=mut.per.yr)) +
  geom_ribbon(fill="lightblue", color="lightblue") +
  geom_errorbar(data=post2, aes(ymin=lower, ymax=upper, x=x),
                inherit.aes=FALSE) +
  geom_point(data=post2, aes(x=x, y=mean), inherit.aes=FALSE,
             shape=21, fill="white") +
  geom_line() +
  facet_wrap(~case) +
  coord_flip() +
  theme(panel.background=element_rect(fill="white", color="black"),
        axis.text=element_text(size=11),
        axis.title=element_text(size=13),
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=13)) +
  scale_x_continuous(breaks=c(1, 3, 5, 7, 9),
                     labels=c(1, 3, 5, 7, 9),
                     limits=c(1, 12)) +
  xlab("Mutations/year\n") +
  ylab("Number of years for IM to GC progression") +
  ylim(c(0, 17))

out_name <- paste0( out_path , "/timing_molecular_clock.pdf" )
ggsave( out_name , p , width=10, height=8)

##########################################################################################
## Average of patient-specific means
average <- round(mean(post2$mean), 1)
